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Abstract 

In this paper, we study the accuracy of adaptively chosen, mapped polynomial ap- 
proximations for functions with steep gradients or discontinuities. We show that, for 
steep gradient functions, one can obtain spectral accuracy in the ongina coor ina ® 
system by using polynomial approximations in a transformed coordinate system with 
substantially fewer collocation points than are necessary using polynomial expansion 
directly in the original, physical, coordinate system. We also show that one can avoid 
the usual Gibbs oscillation, associated with steep gradient solutions of hyperbolic pde s, 
by approximation in suitably chosen coordinate systems. Continuous, high gradient, 
solutions are computed with spectral accuracy (as measured in the physical coordinate 
system). Discontinuous solutions associated with nonlinear hyperbolic equations can 
be accurately computed by using an artifici al viscosity chosen to smooth out the solu- 
tion in the mapped, computational, domain . Thus, we can effectively resolve shocks on 
a scale that is subgrid to the resolution available with collocation only m the physical 
domain. Examples with Fourier and Chebyshev collocation are given. 
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1 Introduction 

During the past decade, spectral methods have become the numerical method of choice 
for the approximate solution of time dependent partial differential equations in many 
branches of computational science. This is due to their ability to yield highly accurate 
approximations of smooth solutions with substantially fewer grid points than would be 
required by comparative finite difference methods. Thus, for complex flows requiring 
detailed simulation, the success of spectral methods has been striking. This is particularly 
true for simulation of turbulent flows and the transition to turbulence (see references in 
[8]). Spectral methods have also become the primary numerical method in global weather 

modeling [14,18]. 

Straight forward application of the spectra! method to problems with highly localized 
phenomena that are not well resolved by the global discretization produce incorrect, glob- 
ally oscillating, solutions due to the Gibbs phenomenon [8,11,27]. This phenomenon shows 
up in linear and nonlinear equations and includes both computational discontinuities (con- 
tinuous solutions with steep gradients that axe not resolved by the finite discretization) 
and real discontinuities (such as shock waves). Various methods [1,8,9,10,12,16,17,20,23] 
have been proposed to extract the correct physical solution from the unwanted oscillatory 
one. Most of these add a finite order truncation error. A spectrally accurate physical space 
filter [1,9,12] and a spectrally accurate vanishing viscosity [22,26] have been demonstrated. 
Both of these approaches, however, perform poorly when used with a realistic number 
of collocation points. More sophisticated, nonoscillatory spectral methods are currently 
under investigation [7]. 

Adaptive grid methods have been used successfully with finite difference and finite ele- 
ment methods for many problems with highly localized phenomenon (see references in [2]). 
The dynamic placement of points to increase local resolution in adaptive finite difference 
and finite element methods is not applicable to spectral methods since the location of 
the collocation points for pseudospectral methods is intimately connected with the global 
convergence properties of the polynomial approximation. (The accuracy of polynomial 
expansions for various choices of collocation points is investigated in [25].) 

An alternative approach to achieving high resolution locally is to use coordinate trans- 
formations which stretch out the region of high gradients so that, in the new coordinate 
system, the high gradient phenomena is well resolved. One can then use standard spectral 
methods to solve the mapped equation in the new coordinate system. Such a technique has 
been successfully applied to a system of reaction diffusion equations arising in combustion 

in [4,5] and independently in [13]. 
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By using adaptively chosen coordinate transformations and computing in the compu- 
tational domain, we are no longer using polynomial approximation for the solution in the 
physical coordinate system. In this paper, we examine the accuracy of mapped polynomial 
approximations. We shall see that, for steep gradient functions, one can obtain spectral 
accuracy in the physical space by using polynomial approximation in a transformed co- 
ordinate system with substantially fewer collocation points than would be necessary with 
polynomial approximation directly in the physical space. We also show that one can avoid 
the usual Gibbs oscillation associated with discontinuous solutions of hyperbolic pde’s, by 
approximation in suitably chosen coordinate systems. Continuous, high gradient solutions 
are computed with spectral accuracy (as measured in the physical coordinate system). 
Discontinuous solutions associated with nonlinear hyperbolic equations can be accurately 
computed by using an artificial viscosity chosen to smooth out the solution in the com- 
putational domain where the solution is already more smoothly represented than in the 
physical domain. Thus, we effectively resolve the shock on a scale that is subgrid to the 
resolution available with collocation only in the physical domain. 

In the next section, we describe the adaptive spectral method and present numerical 
experiments in section 3. Conclusions are presented in section 4. 

2 Numerical Method 

Review of Standard Pseudospectral Method 

We illustrate the standard pseudospectral method [8,11,27] by considering the one dimen- 
sional equation 

u t = G(u) (i) 

where G is a differential operator on the interval x € [a, 6]. The function u(x,t) is approx- 
imated (suppressing the time dependence) by a polynomial 

N 

P N u(x) (2) 

y=o 

that interpolates u(x) at N+l distinct points z 0 , . . . x N at time t. At the collocation points 
x k , if i>j(x k ) = 6 jk , then a y = u(z y ). The expansion polynomial then has the form 

N 

u n{x) = P N u{x ) = ^2 (3) 


2 


W here ^(x.) = V Differentiation of « ia accomplished by analytically differentiating the 

interpolating polynomial ^ 

du N /dx = £ u(xy)0'(x) ^ 

Hence, given the values of u(x) at the collocation points s 4 , we can ffnd the approximate 
derivative (at the collocation points) by a matrix D multiplying the grid function «(*.)■ 
D can be computed explicitly from D = W, = *,(*.)■ The cost of computing the 
derivative is 0(N*) operations from the matrix multiplication. This can be very efficient 
on a vector or parallel computer. Higher derivative matrices can easily be constructed 
from the relation D t = (D)‘. Boundary conditions can also be directly incorpora e in o 

the differentiation matrices. . , fo ., 

Common choices for collocation points and interpolating polynomials [8,27] are 

fourier collocation points 

2irj 


JV + 1 


j = N{Nodd) 


with 


1 . ,( W + !)(*-*,•), .., *-*1 1 

M*) = iv + x "" 1 2 11 2 1 


for periodic problems in x £ [0,2x1, and the Chebyshev collocation points 

Xj = cos nj/N j = 0, . . . , N 


(5) 


with 


T n ( x) = cos(n cos x (x)) 


and 


0j(*) = 


(-lV +1 (l-x 2 )r^(x) 


( 6 ) 


Cj N*{x - Xj) 

where c„ = c N = 2, and c y = l,for 1 < j < N - 1 for nonperiodic problems in x€ [-Ml- 
These particular choices enable the interpolating polynomials to also be written as an 
expansion in orthogonal basis functions {«.)> such that the sums (3 and (4) can be 
computed by an FFT. This reduces the cost of computing the spectral derivative 
O(JV’) to 0(N logJV) operations. The matrix multiplication, however, has the advantage 

that it is more flexible and easier to vectorize. 

Thus, the pseudospectral approximate solution «y = «(*;, t) at time t- is o a, 
by first evaluating the differential operator G at the collocation points at time t-t 
either an FFT or a matrix multiplication, and ;hen updating (in time ) the resulting system 

of ordinary differential equations. 
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Description of Adaptive Procedure 

We now describe the adaptive procedure (first introduced in [5] for a system of reaction 
diffusion equations modeling solid fuel combustion). 

A family of transformations 

5 = g(x, a) ^ 

mapping the interval [a,b] onto itself, and parameterized by the vector o, is introduced. 
Each value of a defines a new coordinate system to which equation (1) can be transformed 
an then solved by the standard pseudospectral method ( in the new coordinate system) 
n he experiments descnbed below we use coordinate transformations with two degrees 
of freedom. One parameter 0 controls the magnitude of the coordinate stretching (or 
compression) about a point determined by the second parameter . For periodic problems 
we map the physical z coordinate system (x € (0,2,r|) to the computational s coordinate 
system (s e [0,2*-]) with the fractional linear transformation 


e” = 


1 — /?e'( 2- 'd 


( 8 ) 


where \0\ < 1 and 0 < 7 < 2ir . For nonperiodic problems, we map the physical coor- 
tr2Hl _M1, the C ° mPUtatiraal " l-MI. with the composite 

s = 2±±)]]. (#) 


* 0 ” ‘2 v 7 x + . 

Tnd different H < '' b ° th “ analytk and •“* i-erted 

Differentiation matrices are easily constructed for the collocation points in the compu- 
tational domain since au/ax = = ,'B/Bs as QD = D, where Q is a diagonal matrix 

whose elements contain the derivative of the map, and D is the differentiation matrix for 
the collocation points in the physical domain. Higher derivatives are obtained by taking 
powers of D. The transformed equations can be easily generated with the matrix formula- 
tion, makmg the exact representation in a particular coordinate system transparent. We 
simp y interpolate u(x,t) to u(s,t) and use D instead of D in (l). 

A particular coordinate system for a function u(x) can be found by minimizing a 

two parameter functional bounding the interpolation error in the maximum norm The 
pointwise interpolation error 


|u(x) -u^x)! = | £ tt^.( z )| 
j=N 


( 10 ) 
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can be reduced by finding a coordinate system n which the leading term in the interpo- 
lation error is minimized. In the nonperiodic case, this leads to the functional derived in 


[4] 


m = 


ds 


y/l — s 2 


In the periodic case, a similar derivation (see Appendix), leads to the H 1 norm 

1(a) = f [H 2 + |u| 2 ]ds. 

Jo 

We will also consider the effects of using the Hj norm, 


( 11 ) 


( 12 ) 


I(ci) = / 2 V|uJ 2 + B \ u *\ 2 + C\u\ 2 ]ds. 
Jo 


( 13 ) 


The algorithm is therefore: 

i) discretize the right hand side of (1) and integrate the resulting ODE’s to t = t x in an 
initial coordinate system. 

ii) find a min which minimizes the monitor function 1(a). 

iii) Interpolate the solution to the collocation points in the new coordinate system using 
the interpolation formula (3). 

iv) Integrate to t = t 2 . 

v) go to (ii). 


3 Numerical Experiments 

In this section we consider the efficiency and accuracy of mapped polynomial expansions. 
We study applications to approximation theory and hyperbolic partial differential equa- 
tions. Using the previously described adaptive procedure, a function is approximated by 
orthogonal polynomials in some transformed coordinate system ’s’. However, we measure 
the accuracy by computing the error in the original ’physical’ coordinate system. 

Approximation Theory 

We measure the approximation error by considering the pointwise error at a sequence of 
points. Since the error at the collocation points is zero by construction, we interpolate to 
a sequence of N' uniformly spaced points, y The pointwise error is given by 

e(yj) = \fiVj) ~ fN{Vj)\ j = 0, . . . , N' ( 14 ) 

where f N (x) is a polynomial of degree N. Next we minimize the monitor function and 
determine a coordinate system, s. We then represent f(s) at the collocation points in the 
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new coordinate system and interpolate to the points s, = s(y y ), the image of the points y y . 
The pointwise error is again measured in the original ’x’ coordinate system. 

In figure la, we display the log of the maximum pointwise error for f(x) = tanh(Mx) 
with M = 8 , 16, 32 in the interval x G [—1,1]. We approximate f(x) by Chebyshev 
polynomials of degree N, where N/M = 1,..,4. The top three graphs are the result of 
standard Chebyshev polynomial approximation, while the bottom three are the result of 
Chebyshev polynomial approximation in a coordinate system selected by the minimization 
of the monitor function ( 11 ). Although both cases converge exponentially, we see the 
mapped polynomials require considerably fewer points to achieve very high accuracy. Since 
the high gradient region occurs in the center of the domain, simply adding Chebyshev 
collocation points provides increasing resolution faster at the boundaries than at the center. 
Therefore it is more effective to use the coordinate stretching to increase local resolution 
than to simply add more collocation points in the original coordinate system. The adaptive 
procedure provides a coordinate stretching where the gradient is largest and maps it to a 
coordinate system where a significant proportion of the collocation points contribute to its 
resolution. 

Next, we examine the convergence properties of the adaptive method when the steep 
gradient region varies away from the center of the domain, where the resolution is coarsest, 
to the boundary, where the resolution is finest. In figure lb. we show the log of the 
maximum pointwise error for f(x) = tanh[M(x - x 0 )] with M = 8,32,128 and x 0 varying 
between the center (x 0 = 0 ) and the boundary (x 0 = 1 ). In this case we use Chebyshev 
expansion with 32 modes. For a given value of M, the top plot is the standard Chebyshev 
polynomial and the bottom plot is the mapped polynomial interpolation. We see that for M 
= 8 and 32 the error is several orders of magnitude smaller with the mapped polynomial 
expansions than with the standard Chebyshev polynomial interpolation. However for 
M— 128, there is insufficient resolution to resolve the gradient even with the mapping; and 
hence the errors are only slightly better with the mapping. The oscillation in the adaptive 
cases is due to the fact that in the discrete search of parameter space (/ 3 , 7 ) to find a 
minimum for ( 11 ) we only looked at values of 7 that were multiples of . 1 , whereas the 
point of maximum gradient (x 0 ) was allowed to vary in multiples of .05. Thus, when x 0 
occurs between multiples of . 1 , the optimal parameters were not found. Still the results 
are much better than if no mapping had taken place. 

Partial Differential Equations 

In the previous section, we have shown that for functions with large gradients, mapped 
polynomials can yield several orders of magnitude more accuracy than the standard spec- 
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tral expansions. We now consider the accuracy of mapped polynomial expansions in the 
context of hyperbolic partial differential equations where the coordinate transformations 
can vary dynamically in time. We consider several prototype examples of wave propa- 
gation including a linear wave equation for advection of a sharp spike and the nonlinear 
Burgers equation for shock formation and propagation. In both cases we will consider 
periodic boundary conditions and hence we use triginometric polynomials for the spatial 
discretization. In order to control the temporal errors in all of our examples we use the 
LSODE solver from the ODEPACK package [6,15] to integrate the resulting system of 
(spatially discretized) ordinary differential equations. This solver uses dynamically chosen 
time stepping in order to guarantee a prescribed error tolerance at each output time. The 
error tolerance was set to 10 -7 . 

In practice we do not change coordinate systems every time step. Instead, we evaluate 
the monitor function /(a) at each time step and compare it with its previous value when 
we last changed coordinate systems. We adapt to a new coordinate system when the ratio 
falls below .7 or above 1.3 . 

Linear Hyperbolic Equations 

Consider the one dimensional linear advection equation 

u t + u x = 0 (15) 

with periodic boundary conditions on the interval 0 < x < 2ir. Nonsmooth initial data 
u(x,0)=f(x) yields the nonsmooth solution u(x,t)=f(x-t) (extended periodically) for all 
t > 0. We choose a sharply peaked Gaussian initial function 

f{ X ) = e -(*-') a /4A* a (16) 

This yields a Gaussian with width 2 Ax, where Ax = 2n/N. If we choose 4 Ax 2 = .025 
then Ax = .078, corresponding to N = 2i r/ Ax « 80 points. Thus we need a local 
resolution corresponding to N=80 points near the peak. In our example we will use 32 
points, therefore the function is not adequately resolved near the peak and behaves as 
if there is a discontinuity. We expect to see global oscillations. In figure 2a, we display 
the computed solution at time t= 0, 15, 30, 45, 60, 64 At. In these examples we use a 
time step of At = .5Ax, and integrate in time with an explicit Adams scheme from the 
LSODE package (MF=10). Note that the LSODE package uses an internally chosen At 
to guarantee the prescribed accuracy. Thus At is only used for output purposes. We 
see that when the solution is not centered at a grid point, the spectral approximation 
displays high frequency oscillations. This could have a serious effect in more complicated 
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systems where the advected quantity is then fed into other components of the system. We 
also point out the observation that the information about the correct physical solution is 
still contained in the oscillatory solution; and the spectral method is able to accurately 
reproduce the solution in future time, when the spike is centered at a grid point. The 
explanation for this is due to the fact that since the interpolation error of the initial 
function is constructed to vanish at the collocation points, and since the spectral method 
introduces almost no dissipation, whenever the solution is merely a translate of the initial 
function to a center about another collocation point it displays no interpolation error. 
Hence we can reproduce, without oscillation, the exact solution at the collocation points 
at times when At is a integer multiple of Ax. At other times when the solution is a 
translate of the initial function between collocation points, we observe the expected global 
oscillation. 

In the adaptive case, we begin by chosing an appropriate coordinate system to ade- 
quately represent the initial function. We therefore minimize the H 1 norm of the initial 
function (16), which we assume to have with as much resolution as needed. The initial 
function is then prescribed at the N fourier collocation points in the new coordinate system. 
The transformed equation 

Ut ■+- (x t )u t = 0 (17) 

is then updated by the standard fourier collocation method. This procedure is then re- 
peated every time step, with the adaption occurring as described previously. In figure 
2b, we display the results of the computation, based on adaption with the H 1 norm, at 
the same output times as in the non-adaptive case (figure 2a.). For output purposes, we 
have interpolated the solution back to the physical x coordinate system. In figure 2c, the 
same computation was done using the H 2 norm of u as the minimizing functional. We also 
present, in figure 2d, for comparison, results of the adaptive procedure with the spatial 
differentiation matrix replaced with the fourth order finite difference matrix. 

We compare these cases by computing the I 2 errors of the grid functions ie 


HOII2 


1 N 

{ xr ) ^ |^eioct(^j> f) ^'computedi^'j ^ 
1 3 = 0 


and the /<*, error 

li e (OI|oo = Tnax j\ u exact{ x jj t) — U comptl ted{Xj , t) |. 

The maximum time errors 


max 
< e [ o,2k] 


KOII 


2,00 


for the four schemes is presented in table 1. 


(18) 


( 19 ) 
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Method 

max I 2 error 

max loo error 

Standard fourier collocation 

.038 

.096 

Adaptive fourier collocation using H 2 norm 

.003 

.009 

Adaptive fourier collocation using H 1 norm 

.006 

.012 

Adaptive 4-th order finite diff. using H 2 norm 

.143 

.446 


These results suggests that adaption based on minimizing the H 2 norm is preferred 
over adaption based on the H * norm. ^Vhen the sole criteria for choosing the coordinate 
system is the H 1 norm too much weight is placed on the steep gradient (high frequency) 
region of the solution at the expense of the smooth (low frequency) region. We can see in 
figure 2b, the presence of a low frequency error component in the smooth, constant, portion 
of the solution. The use of the H 2 norm provides a smoother coordinate mapping with 
more emphasis on the smooth region of the solution. Thus, the low frequency components 
of the solution is more accurate. Both adapted solutions, however, provide significant 
improvement over the nonadapted solution and the adapted fourth order finite difference 
solution. 

Nonlinear Hyperbolic Equations 

As a model for the formation of discontinuous solutions, we consider Burgers equation 

u, + i(u 5 ), = 0 (20) 

with periodic boundary conditions on the interval 0 < x < 2 ir. Smooth initial data 
u(x,0)=f(x) will develop discontinuous (weak) solutions in finite time. An entropy condition 
is required to pick out the correct physical solution [21]. 

The fourier collocation method will break down and yield global oscillations as soon as 
the solution steepens up too much to be resolved on the finite grid. The adaptive procedure 
will delay this breakdown, but its occurrence is inevitable. Artificial viscosity can smooth 
out the solution so that it becomes resolvable on the finite grid. This introduces a finite 
order truncation error and therefore we loose the spectral accuracy. Recently [22,26], a 
spectrally accurate vanishing viscosity solution has been shown to converge to the entropy 
solution of (20). However, computations with realistic number of collocation points still 
shows global oscillations [26]. 

We propose here to couple the adaptive pseud ospectral method with the artificial vis- 
cosity method to produce oscillation free solutions using significantly less artificial viscosity 
than would be necessary with standard fourier collocation. This will reduce the truncation 
error introduced by the artificial viscosity. Gottlieb et al. [10] have shown that spectral 
collocation methods are capable of resolving a shock over one effective grid interval. We 
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therefore expect, that by appropriate coordinate stretching, the adaptive pseudospectral 
method can resolve shocks over a significantly smaller effective grid interval; the equivalent 
of using several times more collocation points. 

We approximate the solution of the inviscid Burgers equation (20) by the viscous version 


u t + u 2 ) x = uu x 


( 21 ) 


where uu xx is the artificial viscosity term, and v is a grid dependent viscosity coefficient. It 
is well known [21], that in the limit u J. 0 solutions of (21) converge to the entropy solution 
of (20). In our numerical experiments we use initial data f(x) = sin(x). The inviscid 
Burgers equation develops a shock discontinuity at x = 7 r. Again, the fourier collocation 
differentiation matrix is used for the spatial differentiation. The resulting system of ODE’s 
was integrated using an explicit Adams method (MF=10) from the LSODE package. The 
solution was integrated up to time t = 1.47, after the shock has formed. We determined 
experimentally that for collocation with 32 points, a viscosity oi v = .08 is needed to 
obtain an oscillation free solution. The adaptive solution based on the H 1 norm only needs 
a viscosity coefficient of v = .02 to yield an oscillation free solution. In figures 3a,b,c we 
show a) a time history of the solution up to t=1.47, b) the exact and computed solutions at 
and c) a time history of the errors for the non adaptive fourier collocation solution 
of (21) with artificial viscosity coefficients of u = .08. In figure 4a,b,c we show the same 
plots for the adapted solution with an artificial viscosity coefficient v = .02. 

Comparing the two cases we see the non adapted solution requires four times as much 
viscosity to reduce the global oscillations. The effect of this larger viscosity can be seen in 
a much larger region about the shock and, in fact, prevents the computed solution from 
shocking up. The maximum error occurs just before the solution shocks up, when the 
second derivative is greatest. The smaller viscosity used in the adaptive solution accounts 
for a reduction in the maximum pointwise error by a factor of four. The larger viscosity 
needed for the nonadaptive solution prevents the solution from steepening any further, 
and the maximum pointwise error persists after the shock forms. However the adapted 
solution, with the smaller viscosity, is able to steepen further and approach the true shock 
solution. After the shock forms, the solution is essentially piecewise linear and the effect 
of the second derivative truncation error is reduced. The error is then reduced another 
factor of four. 

Finally, we consider the case of a moving shock. If we use the initial condition f(x) = 
.3 + .7sin(x), the solution shocks up and starts starts to propagate. In figure 5a, we 
plot a sampling of the time history of the solution up to t=4.1. The pointwise error at 
these sampled times are plotted in figure 5b. As in the previous case, we use adaption 
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based on the H 1 norm and a viscosity coefficient of u = .02 . We see that, as before, the 
maximum pointwise errors occurrs just before the solution shocks up, when its curvature 
is maximum. After the solution shocks up it almost piecewise linear and the truncation 
error due to viscosity is reduced. The solution then propagates without exhibiting global 
oscillations. Subsequent peaks in the errors (displayed in 5b) is due to the discrete search 
of parameter space in the minimization of the monitor function (12). More sophisticated 
minimization routines will be used in the future. 


4 Conclusions 

We have demonstrated that mapped polynomial approximation can provide significant 
improvements in accuracy and efficiency for approximating functions with steep gradients. 

Using analytic, dynamically varying, simple coordinate transformations we achieve im- 
provements in convergence of several orders of magnitude for continuous functions with 
steep gradients. These results are robust in the sense that even with crude minimization 
procedures, that don’t find the optimal coordinate transformation, we achieve significant 
improvements over the nonadapted case. 

Discontinuous functions require smoothing in addition to coordinate transformations. 
In this work we have used artificial viscosity to provide the smoothing. Although this 
adds a finite order truncation error, the effects are greatly reduced by computing in the 
transformed coordinate system where less artificial viscosity is needed. Artificial viscosity 
does have the drawback that it reduces the allowable time step for explicit methods from 
0(N 2 ) to 0(JV 4 ) for Chebyshev based collocation. We are currently working with alternate 
smoothers based on filtering the expansion coefficients in spectral space [10,23]. These 
results will be reported elsewhere. 

In this paper, we have restricted our attention to problems with only one high gra- 
dient region. Problems with multiple steep gradient regions can be treated using more 
complicated mapping functions. This will increase the dimensionality of parameter space 
searched in the minimization of the monitor function. A better approach is to couple the 
adaptive spectral method with a multidomain approach [19,24] so that the simple map- 
pings described above can be used in each domain. We are currently investigating this 

approach. 
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Appendix A 

We now derive a bound for the maximum norm of the interpolation error of a periodic func- 
tion. Using the standard notation, we represen t a periodic function u(x) as an expansion 
of trigino metric basis functions 

«(x) = f) a,e“- (22) 

Jk=-oo 


where 


a* = — [ \(x)e~ ikx dx, |fc| = 0, 1, 2, 
2 7T JO 

If we use the the galerkin approximation of u(x). 


N 


U N (x) = Pnu{x) = Y a * e 

k=-N 


ikx 


(23) 


(24) 


then the pointwise error is 


u 


(x) — Ujv_l(x)| I Y/ afce I* 

|fc|>N 


We can now bound the spectral interpolation error in the maximum norm by 


max 0 < z <2r\u{x) — ujv_i(x)| < Y l°*l* 

\k\>N 


(25) 


(26) 


Applying integration by parts to (23) p times and noting the periodicity of the boundary 
conditions, we have 

_ 1 (-tr 1 [*’ p = 1,2, ... (27) 

2tt (iN) p Jo 


O.N 


We can bound a N by applying the Cauchy-Schwartz inequality to (27) 

|a»|<i{/ o 2 V>| ! <te>’ P=1.2,-.. 


(28) 


It is therefore reasonable to use the functioned 

1(a) = A|u„| 2 + B|u,| 2 + C\u\ 2 ds (29) 

as the measure for the interpolation error in a coordinate system s. 

Although this is the same functional used in [5,13] for Chebyshev polynomial approxi- 
mation, the estimates leading to (29) in that case are not sharp. In fact, (29) was reported 
in [4] to not be sufficiently sensitive to variations in the behavior of the solution. When 
dealing with periodic functions, however, the same type of sharp estimates that lead to 
the alternate functional used in [4], here lead to the H* norm functionals (eqs. (12) and 

(13)), as used originally in [5,13]. 


13 


References 


[1] Abarbanel, S., Gottlieb, D., and Tadmor, E.; Spectral methods for discontinuous prob- 
lems in Numerical Methods for Fluid Dynamics, ed. by K. W. Morton and M. J. 
Baines, Oxford University Press (1986), pp. 129-153. 

[2] Babuska, I., Chandra, J., and Flaherty, J.; Adaptive Computational Methods for Par- 
tial Differential equation , SIAM, (1984). 

[3] Basdevant, C., Deville, M., Haldenwang, P., Lacroix, J., Orlandi, D., Quazzani, J., 
Patera, A., and Peyret, R.; Spectral and finite difference solutions of the Burgers 
equation , Computers and Fluids 14 (1986), pp. 23-41. 

[4] Bayliss, A., Gottlieb, D., Matkowsky, B., and Minkoff, M.; An adaptive pseudospectral 
method for reaction diffusion problems, ICASE report No. 87-67 (1987), 40 pages. 

[5] Bayliss, A. and Matkowsky, B.; Fronts, relaxation oscillations and period doubling in 
solid fuel combustion, J. Comp. Phys. 71 (1987), pp. 147-168. 

[6] Byrne, G. D. and Hindmarsh, A. C.; Stiff ODE solvers: A review of current and 
coming attraction, J. Comp. Phys. 70 (1987), pp. 1-62. 

[7] Cai, W., Gottlieb, D., and Shu, C. W.; Non-oscillatory spectral fourier methods for 
shock wave calculations, ICASE report No. 88-37 (1988), 31 pages. 

[8] Canuto, C., Hussaini, M. Y., Quarteroni, A., and Zang, T. A.; Spectral Methods in 
Fluid Dynamics, Springer- Verlag (1988), 512 pages. 

[9] Gottlieb, D.; Spectral methods for compressible flow problems in Proc. 9th Int. Conf. 
Numerical Methods in Fluid Dynamics, ed. by Soubbarameyer, J. P. Boujot, Springer- 
Verlag (1985), pp. 48-61. 

[10] Gottlieb, D., Lustman, L., and Orszag, S.; Spectral calculations of one-dimensional 
inviscid compressible flows, SIAM J. Sci. Stat. Comput. 2 (1981), pp. 296-310. 

[11] Gottlieb, D. and Orszag, S.; Numerical Analysis of Spectral Methods: Theory and 
Applications, SIAM (1977). 

[12] Gottlieb, D. and Tadmor, E.; Recovering pointwise values of discontinuous data within 
spectral accuracy, in Progress and Supercomputing in Computational Fluid Dynamics, 
ed. by E. M. Murman and S. Abarbanel, Birkhauser (1985), pp. 357-375. 


14 


,r,. 0 nthc use of spectral methods for the numerical solution 

1131 “iTms cZ, AP, — ■ »* - («*>> - ^ 

N a , - w— ■ - - “ ^ “ 

John Wiley and Sons (1980). 

v a n-ODEFACK a systematized collection of ODE solvers, 

1151 ^ 
m m : y -’ k °? ^ U985), 

Euler equations, part 1: Courier me 

PP 64 7 „ „ • n Salas M D.. and Zang, T. A.; Spectral methods for the 

l»l U t : Y ; yj cli sTeu methods and shade fitting AIAA J. 23 (1985), PP- 

Euler equations, part d. oneo yon 

234-240. 

. A p M .The use of spectral techniques in numerical weather 
1181 Computations in Fluid Mechanics, Lectures in Applied 

Mathematics 22 (1985), pp- 1-41- 

„ " “ — " 

Appl. Numer. Math. 2 (1986), pp. 221 241. 

11 (1972). 

. _ . E . Ana lysis cl the spectral vanishing viscosrty method for 

[22) Maday, V, and Tadm r, E '' A ^ ’ 88 . 4 (1988) , 34 pages. 

periodic conservation laws, ICASL Keporr 

, . , 0sher s.. The fourier method for nonsmooth emteal 

[23] Majda, A., McDonough, a ' 

data , Math. Comp. 32 (1978), PP- 1041-1081. 

*„ d Street C L • Improvements in spectral collocation through a mul- 

[24l Macaraeg, M., and Stre , • (1986) pp. 95-108. 

tiplc domain technique, Appl. Numer. Math. 2 (1986), pp 

_ . nd Turkel E • Global collocation methods for approximation and the 
[251 Solomonoff, A. and Turkel, t , gg _ 60 86)j 65 pages . 

solution of partial differential eguatiom, ICASE Repo 


15 



1261 

ential Equations, SIAM ( 198 ^ ITp^es. J" ^ '°' ” afft r ‘ 


16 



Figure Captions 

Figure la. Log of pointwise error for fn{x) « tanh{Mx) where jj = 1,2, 3, 4. 


Figure lb. 


Figure 2a. 


Figure 2b. 
Figure 2c. 


Log of maximum pointwise error for fsz^x) « tanh[M{x £<})]• M 
8,32,128 and 0. < x 0 < 1. 

Computed solution for the linear advection equation (15) with 32 collo- 
cation points at time t = 0,15,30,45,60 and 64 At, with At = .5(|f). 

Same as (2a) except solution adaptively computed based on H 1 norm. 
Same as (2a) except solution adaptively computed based on H 2 norm. 


Figure 2d. 


Figure 3a. 


Figure 3b. 


Same as (2a) except solution adaptively computed based on H 2 norm 
and fourth order spatial differentiation is used. 

Time history for viscous nonadapted Burger equation up to t=1.47. Vis- 
cosity coefficient is u = .08. 

Computed and exact solution at time t = 1.47. 


Figure 3c. Pointwise error history for solution displayed in 3a. 

Figure 4a. Time history for adaptively computed viscous Burger equation up to t = 

1.47. Viscosity coefficient is u = .02. 

Figure 4b. Computed and exact solution at time t = 1.47. 


Figure 4c. Pointwise error history for solution displayed in 4a. 

Figure 5a. Sampling of time history for adaptively computed viscous Burger equa- 

tion up to t = 4.1. Viscosity coefficient is v = .02. Initial function 
u(x) = .3+ .7stn(x) 

Figure 5b. Pointwise error history for (solution displayed in 5a. 
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